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ABSTRACT 



We study the dynamical evolution of globular clusters containing primordial 
binaries, including full single and binary stellar evolution using our Monte Carlo 
cluster evolution code updated with an adaptation of the sing le and binary stel- 
lar evolution codes SSE/BSE from iHurley et al.l (120001 . 120021 ). We describe the 
modifications we have made to the code. We present several test calculations and 
comparisons with existing studies to illustrate the validity of the code. We show 
that our code finds very good agreement with direct iV-body simulations includ- 
ing primordial binaries and stellar evolution. We find significant differences in 
the evolution of the global properties of the simulated clusters using stellar evo- 
lution compared to simulations without any stellar evolution. In particular, we 
find that the mass loss from stellar evolution acts as a significant energy produc- 
tion channel simply by reducing the total gravitational binding energy and can 
significantly prolong the initial core contraction phase before reaching the binary- 
burning quasi steady state of the cluster evolution as noticed in Paper IV. We 
simulate a large grid of clusters varying the initial cluster mass, binary fraction, 
and concentration and compare properties of the simulated clusters with those 
of the observed Galactic globular clusters (GGCs). We find that our simulated 
cluster properties agree well with the observed GGC properties. We explore in 
some detail qualitatively different clusters in different phases of their evolution, 
and construct synthetic Hertzprung-Russell diagrams for these clusters. 



1. Introduction 



Star clusters in general, and the Galactic globular clusters in particular, have been 
studied extensively for many years. As tracers of the Galaxy's potential, their dynamical 
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history tells us something about the formation and evolution of our Galaxy. As dense stellar 
environments, their interesting constituent populations (including, e.g., blue straggler stars, 
cataclysmic variables, and low-mass X-ray binaries) inform our understanding of binary stel- 
lar evolution through its interaction with dynamics. The study of the evolut ion of globular 



and o ther dense star clusters has had a somewhat long and varied history (e.g.. lHeggie fc Hut 



20031 ). Before observations showed that globular clusters contained dynamically significant 
numbers of binaries, theoretical efforts focused on understanding the process of core collapse 
and the ensuing post-collapse evolution driven by three-body binary formation. Once it 
became clear in the early 90s from observations that clusters contained sufficient numbers of 
binaries that they must have been born with substantial "primordial" populations, theory 
emphasized properties of clusters in the "binary burning" phase, in which the cluster core 
is supported against collapse by super-elastic dynamical scattering interactions of binary 
stars. More recently it has been realized that pure point-mass interactions of binaries re- 
sult in equilibrium cluster core radii in the binary burning phase that are a factor of ~ 10 
smaller than what is observed, and so many efforts have focused on alternative cluster en- 
ergy sources such as central intermediate-mass black holes, expedited stellar mass loss from 
compact object formation via collisions, or prolonged mass segregation of compact objects 



(IHeggie et al.ll200fil : iFregeau fc Rasioll2007l : lTrentill2006l : IChatterjee et al.ll2008l : Merritt et al. 



2004 ). In a similar vein, recent theoretical work, combined with observations showing that 
the core binary fraction in many clusters is fairly low (< 10%), s uggest that clusters m ay be 
born with remarkably low binary fractions of just a few percent (IFregeau et al.ll2009l ). Such 
a small primordial binary fraction would be surprising since observations of young stars sug- 
gest that star formation yields binary fractions of order ~ 50%. Clearly, our understanding 
of globular cluster evolution has changed considerably over the past few decades, and much 
of that change has been driven by numerical simulations. 

Among computational tools for st udying the dy namical evolution of star clusters, the 
Henon Monte Carlo (MC) technique (Henon 1971al lbl) represents a balanced compromise 
between exactness and speed. The MC method allows for a star-by-star realization of the 
cluster, with its N mass shells representing the N stellar objects in the cluster (either single 
or binary stars). It assumes, most importantly, spherical symmetry and diffusive two-body 
relaxation, allowing time integration on a relaxation timescale, and a computational cost that 
scales as N log N. We have developed our Henon MC cluster evolution code (which we call 



CMC for "Cluster Monte Ca rlo") over the past decade (IJoshi et alJl2000l . 12001 



Fregeau et al. 



20031 ; IFregeau fc Rasioll2007l . henceforth Papers I, II, III, and IV, respectively). Since it allows 
for a star-by-star description of the cluster at each timestep, it is relatively easy to add 
physical processes beyond two-body relaxation to the code. We have previously added the 
effects of a Galactic tidal field, dynamical scattering interactions of binary star systems, and 
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physical collisions between stars. In this paper we describe the addition of stellar evolution of 
single and binary star systems. Many stars in a cluster evolve internally on a timescale shorter 
than the age of the cluster. At early times they may lose a substantial fraction of their mass 
via stellar winds. At later times they may evolve off the main sequence, changing their masses 
and radii (and hence collision cross section), and possibly receiving systemic velocity kicks 
when they become compact objects. Since the binding energy of binary stars is an important 
fuel source that can postpone the deep core collapse of star clusters, stellar evolution of 
binary systems directly affects their global evolution. Conversely, the properties of the 
cluster environment feed back on stellar evolution, modifying the evolutionary pathways of 
binary systems and the properties and number s of interesting binary syst ems relative to the 
low-density Galactic field (e.g., X-ray binaries; llvanova et al.ll2006l . 120081 ). 



Previous cluster evolution studies that include stellar evolution have improved our un- 
derstanding of the global evolution of clusters greatly, identifying several distinct stages of 
evolution. At early times, as the stars are forming and the most massive stars have already 
begun nuclear burning, the cluster loses mass through residual gas expulsio n and stellar 
wind s , resu lting in cluster expansion during the first few Myr of evolution Hurlev et al. 



2001 



20021 ). Shortly thereafter, if a runaway collision scenario is avoided (e.g., 



Freitag et al. 



2006al ). two-body relaxation dominates, resulting in a fairly long-lived (from a few to tens of 
Gyr) phase of core contraction. Once the core density becomes high enough for the energy 
generated in binary scattering interactions to balance the energy carried out of the core 
by two-body relaxation, a potentially very long-lived (up to tens of Hubble times or more) 



phase of "binary burning" ensues (e.g.. iHurleyl 120071 ; iFregeau fc Rasiol 120071 ). Once the pop 



ulation of binaries is exhausted in the core, the cluster goes into deep core collapse. In the 
classical, point-mass limit, deep core collapse is arrested by the formation of a " three-body 
binary" and followed by a phase of gravothermal oscillations ( jHeggie fc Hutll2003l ). However, 
three-body binar y formation may be inhibited by stellar collisions in sufficiently young and 
massive clusters ( jFreitag et al.ll2006bl ). 



With the exception of a few recent simulations, most numerical studies that include stel- 
lar evolution have either been limited in the number of stars they can treat or h ave adopted a 
narrow initial mass function with very simplified stellar evolution recipes (e.g., 



Giers 



j pted c 
z 1998 



Giersz fc Spurzemll2000l : Ijoshi et alJl2000L l200ll ; IFregeau et al.ll2003l : IFregeau fc Rasioll2007h . 

Stars in star clusters are born with a range of masses up to ~ 100 M^, and down to at 



least the hydrogen-burning limit (e.g.. ISalpeterlll955l ; iMiller fc Scaldll979i ; lKroupall200ll ). so 
one should evolve the full spectrum of stellar masses as realistically as possible to properly 
treat the influence of stellar evolution on global cluster evolution. Emphasis has recently 
been placed on comparing observed properties of Galactic globular clusters with the oretical 
predictions. Comparison of observed cluster structural properties with theory (e.g., iHurley 



-4- 



20071 ; iFregeau fc Rasioll2007D suggests that either an additional energy sourc e is "puffing" up 
cluster cores (IMackey et al.ll2008l ; iTrenti et al.l 120071 ; IChatterjee et al.l 120081 ), or perhaps th e 
clusters are not in the expected evolutionary states, namely binary-burning (lFregeaull2008l ). 
For example, recent iV-body simulations by iHurleyl (120071 ) show that the core contraction 
phase can last a Hubble time, resulting in a cluster core radius that is larger than one would 
expect were the cluster in the binary burning phase. In these simulations, the core contrac- 
tion phase is prolonged by mass loss from stellar evolution. Clearly, stellar evolution may be 
an important component in globular cluster evolution. 

To more prope rly treat stellar evo l ution , we have recently coupled to CMC the stellar 
evolution recipes of iHurley et al.l (120001 . 120021 . hereafter referred to as BSE). We choose BSE 
for ease of implementation, and for more direct comparisons with iV-body calculations, which 
commonly use BSE for stellar evolution. In §2] we describe the implementation of stellar 
evolution in our code. In £J3] we validate it by comparing with existing cluster evolution 
calculations in the literature. In §4] we demonstrate the importance of stellar evolution by 
comparing simulations that do not include it. In $5] we apply our newly updated code to 
the evolution of a large grid of cluster models, highlighting typical behavior and comparing 
with observations. Finally, in §S]we summarize and conclude. 



2. Method 



CMC treats a number of important physical processes, including two-body relaxation, 
the tidal effects of a host galaxy, strong binary-binary (BB) and binary-single (BS) scattering 
interactions, and direct physical collisions between stars (Paper IV). Here we describe the 
recent addition of BSE (IHurley et al.l l2000l , 120021 ) to treat stellar evolution of single and 
binary stars. 



2.1. Stellar Evolution 



For ease of implementation and for more direct compariso ns with direc t iV-body w e 
use the BSE stellar evolution routines, as described in detail in IHurley et al.l (120001 . 120021 ). 
For single stars, BSE comprises analytic functional fits to theoretically calculated stellar 
evolution tracks as a function of metallicity and mass. Binary star systems are treated using 
the same fitting formulae for each star, but with treatments of physics relevant to binaries, 
including stable and unstable mass transfer, common envelope evolution, magnetic braking, 
tidal coupling, and the effects of neutron star and black hole birth kicks. As described above, 
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CMC allows for a star-by-star description of the cluster at every timestep, allowing for the 
inclusion of additional physics. The stellar properties of binary and single stars are updated 
in step with dynamics via function calls to the BSE library. As described in detail in Papers 
I-IV, CMC uses a shared timestep which must be set as small as the smallest characteristic 
timescale for each physical process. We set the characteristic timescale for stellar evolution 
to the time for the cluster to lose a fraction 0.001 of its total mass. In this way, any cluster 
expansion from stellar evolutionary mass loss is properly resolved. 

One aspect of our BSE implementation requires special mention, however. The evolution 
of high mass stars (> 100 M ) is rather uncertain and can vary greatly depending on the 
wind mass loss prescription. These high mass stars are also quite rare and short lived, so 
observational constraints are limited. In BSE stars with mass > 100 M & are evolved as if 
they were 100 M & stars. When their dynamical properties are returned from BSE, their 
original (> 100 M ) mass is returned. 



2.2. Collisions 



As described in Paper IV, collisions are treated in the sticky-sphere approximation, 
which w as shown to be rema rkably accurate for the velocity dispersions found in globular 
clusters ( Freitag et al. 2006al lbl). When two stars collide, their properties (e.g., stellar type 
an d effective age) are s et according to the BSE merger matrix and prescriptions as described 



m iHurley et al.l (120021 ). In CMC this is implemented in the following simple way. The two 
stars are passed to BSE as a very tight, eccentric binary with nearly zero pericenter distance 
and evolved for a very short time until they merge. The properties of the merger product 
are then naturally set within BSE, and returned to CMC as a single star. BSE by default 
assumes full mixing of nuclear fuel during a collision involving MSSs. We adopt this same 
rejuventation prescription for the simulations presented in this paper, but note that the 
amount of mixing in the collision product should depend on the details of the interaction 
parameters lea ding to the collision, as we l l as the evolutionary s tages of the collision pro- 



genitors (e.g., lLombardi et aI.lll995L 1199ft ISills et al.lll997l . l200a lLombardi et al.ll2002f ). In 



fact, it is found using detailed SPH calculations that the amount of mixi ng as a result of 



a col l ision may be minimal, especially for collisions involving evolved stars (lLombardi et al. 
19951 . Il996h . 
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Fig. 1. — Evolution of the total cluster mass for clusters with various initial numbers of stars 
and masses. In each case a Galactocentric distance tqc — 8.5 kpc and a standard Galactic 
tidal field are assumed. The initial number of stars (N) is noted in each panel. The solid 
black line and the dashed black line in each panel show CMC results with the apocenter 
criterion and energy criterion, respectively. The dotted blue line in the first four panels show 
the NBODY4 results for the same initi al clusters for compari s on. T he NBODY4 data is 
taken from the simulations described in iBaumgardt fc Makinol (120031 ). Similar simulations 
using NBODY4 do not exist for a higher N. In each case the energy criterion loses mass 
at a faster rate than the apocenter criterion. When available, the results using the energy 
criterion agree better with the direct iV-body results. 
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2.3. Tidal Truncation Treatment 



Globular clusters are not isolated systems, but are in fact subject to the tidal field of 
their host galaxy. The assumption of spherical symmetry inherent in MC codes like CMC 
does not allow for a direct calculation of stellar loss at the tear-drop shaped tidal boundary. 
Instead, MC codes employ an effective tidal mass loss criterion that attempts to match 
the tidal mass loss found in direct iV-body simulations. Since stars are lost from the tidal 
boundary on a dynamical timescale and MC codes operate on the (much longer) relaxation 
timescale, the appropriate effective criterion is not obvious. There have been two main 
suggestions in the literature for the appropriate tidal truncation criterion. Perhaps the most 
natural is to immediately strip any star whose apocenter (r a ) lies outside the Roche lobe 
radius of the cluster (which we call the tidal radius, r t ): 



r a > r t . 



X 



This "apocenter criterion" has been used exclusively in CMC previously (Papers I-IV). 
Earlier, in the absence of large- N direct iV-body simulations, comparisons were made with 
2-D Fokker-Plank models and the apocenter criterion showed excellent agreement (for details 
see Papers I and II). 

Another simple criterion is the "energy criterion," in which any star with a specific 
orbital energy above some critical energy is immediately stripped: 



E nr h > l.l 



H ■ 



(2) 



where <j) t is the cluster potential at the tidal radius (jSpitzerl 119871 ). However, a stellar orbit 
that instantaneously satisfies the above criterion may still remain bound if it is scattered 
back to a lower energy orbit before it can escape. To account for this effect a less obvious 
but empirically va lidated correction factor to the energy criterion above is suggested by 
Giersz et all feoosh : 

E OTh > a<p t , (3) 
where a is an iV-dependent parameter given by 



1.5-3 



ln{jN) 
N 



1/4 



(4) 



We have re-examined and tested these two criteria (Equations [T] and [3]) for tidal stellar 
loss to d etermine which one agr e es be tter with the latest results from direct iV-body sim- 
ulations. iBaumgardt fc Makind (120031 ) study in detail the tidal dissolution timescales of a 
cluster in a tidal field varying the initial number of stars and the initial mass of the cluster. 
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We have repeated a large subset of this extensive study of tidal disruption using CMC with 
both tidal truncation criteria, and compared the results. 

We followed th e exact same prescription for setting up the initial conditions of the clus- 



ters as described in lBaumgardt fc Makind (120031 ). Each cluster in this set of runs is assumed 



to have a circular orbit around the Galactic center with radius 8.5 kpc. The Galactocentric 
circular velocity Vq is assumed to be 220km/s. The initial position and velocity of each star 
of the cluster is chosen from a King model distribution function with central concentration 
parameter Wq = 7. The initial masses of the stars are chosen according to the Kroupa 2001 
mass function in the range 0.1 — 15 M Q . We vary the initial number of stars in the clusters 
between 16000 and 10 6 . There are no primordial binaries in these simulations. 

Figure [T] shows the evolution of the bound mass in a cluster in the standard Galactic 
tidal field as described above. We find that the energy criterion results ag ree better with the 



direct iV-body results, when available from iBaumgardt fc Makind (120031 ). The agreement is 



much poorer when using the apocenter criterion. For example, for the direct iV-body run 
with initial N = 64K , the time when the cluster loses 80% of its initial mass is ~ 1 1 Gyr. 
Using CMC with the energy criterion this value is ~ 10 Gyr, whereas, with the apocenter 
criterion the same cluster does not lose 80% of its initial mass within 20 Gyr, the integration 
stopping time. 

Since the energy criterion gives a significantly better agreement with existing direct 
iV-body results (Figure (TJ see also more detailed comparisons in §3]), we adopt the energy 
criterion here, in contrast to what was used in our earlier works (Papers II - IV), unless 
otherwise mentioned. 

One should note, however, that the cluster mass range where direct iV-body results are 
available is not representative of the actual GGC population. Since no direct iV-body results 
exist for more massive clusters, it is not possible at present to determine which approximation 
is more accurate for larger N. 

The two criteria above are treated as initial options in CMC and either one can be 
selected at the beginning of a simulation. At each timestep, the amount of mass lost is 
calculated using the chosen criterion in an iterative way to obtain the bound mass (see 
Paper II for details). 
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3. Comparison with direct iV-body Results 



In this section we validate our treatment of stellar evolution by comparing wit h results 
from previously published studies using the popular direct iV-body code NBODY4 (lAarseth 
20031 ). Since the direct iV-body simulations employ the least degree of assumptions, we tend 
to trust them more for validation purposes. 

One of the biggest simulations treating all rel evant physical effects i ncluding primo rdial 
binaries and stellar evolution was performed by iHurley et al.l (120071 ); iHurleyl (120071 ). In 
particular they studied the evolution of the core properties, binary number fractions in the 
core as well as in the full cluster, and the evolution of the bound number of stars. Since 
both these works present data from a common set of simulations we henceforth collectively 
call them Hurley07. 



3.1. Initial Conditions 

We choose from Hurley07 the simulations with a number iV, = 10 5 initial objects (the 
largest initial iVj in their set of simulations), with primordial binary fractions ff, = 5% and 
fb = 10% (their K100-5 and K100-10 models, respectively). Throughout this work we count 
binary center of mass as one object. Thus a cluster with iV, = 10 5 and = 5% initially 
has 95000 single stars and 5000 binaries. We simulate clusters with exactly the same initial 
conditions using CMC. The initial stellar positions and velocities are chosen fr om a virialized 



Plumm er sphere. The stellar masses are chosen from the IMF presented in iKroupa et al. 



(119911 . henceforth K91) in the range 0.1 - 50 M . Metallicity is fixed at z = 0.001. Each 
cluster has an initial virial radius r v = 8.5 pc which corresponds to a roche-filling cluster 
with tidal radius r t ~ 50 pc, consistent with the standard Galactic tidal field with a Galactic 
rotation speed 220km/s at a Galactocentric distance 8.5 kpc. The companion mass in each 
binary is chosen from a uniform distribution in mass ratio in the range 0.1 — m p M , where, 
m p is the primary mass. The bina ry period is drawn fr om a distribution flat in logarithmic 



intervals in the semi- major axes (jEggleton et al.l Il989l ) and the eccentricities are thermal 



following Hurley07. We call these simulations hcnle5b5 and hcnle5bl0, respectively (Table 



3.2. Comparison of Results 



For the global evolution of a dense cluster the evolution of the core is extremely im- 
portant since throughout the evolution the global properties of the cluster are determined 
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Fig. 2. — Evolution of r c (top) and r c /rh (bottom). The solid black lines show results from 
simulation using CMC and the red dashed lines show results from Hurley07. All data from 
Hurley 07 used for comparison i n this work are extracted using ADS's Dexter data extraction 
applet (IDemleitner et al.ll200ll ). 
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by the balance of energy in and out of the core. The core radius (r c ) is also one of the 
most easily observable structural properties of a cluster. Moreover, the properties and the 
evolution of r c is one of the easiest theoretical way to characterize the distinct phases of a 
cluster's evolution. Thus a basic test for the validity of a cluster simulation is to compare 
the evolution of the core radius (r c ) and the ratio of the core to half mass radius (r c /r^). 

Note that for all our simulated clusters in this work r c is the density-weighted core radius 
(ICasertano fc Hutlll985l ) commonly used in iV-body simulations, unless otherwise specified. 
This i s not a direct ly observable quantity and can differ from the observed r c by a factor of 
a few JiIurlevll2007L 



Figure [2] shows the evolution of r c for run hcnle5b5 (Tabled]) and K100-5 in Hurley07. 
The scale-free quantity r c /r h is also plotted for each run. The core radius expands due 
to stellar evolution mass loss during the first ~ 2 x 10 2 Myr. The core then contracts at 
a steady rate till a little after ~ 1.5 x 10 4 Myr. The core radius then attains a relativel y 
steadier value as the cluster reaches the binary-burning phase (e.g. iFregeau fc Rasiol 120071 ). 
All these qualitatively different phases of the evolution of a cluster is reproduced using CMC 
with excellent agreement. 

One of the key results of Hurley07 is that the overall binary fractio n (fb) remains 



close to the primordial value throughout the evolution of the cluster (also see IFregeau et al. 



20091 ). This result has immense observational significance. In practice only the present day 
properties of a cluster are observed. This result from Hurley07 indicates that if a present 
day binary fraction of the cluster close to can be observed the primordial hard fb should 
have been close to this observed value. Figure [3] shows the evolution of the core (fb,c) and 
the overall fb from CMC simulation hcnle5b5 and direct iV-body simulation presented in 
Hurley07. Binaries preferentially sink to the center due to mass segregation and the single 
stars typically get tidally disrupted from the tidal boundary. These two effects compete 
with each other — the first reduces and the second increases fb outside the core. For the 
simulated cluster these two effects more or less balance each other. For the simulated cluster 



Table 1. Initial conditions for comparison runs with Hurley07 



Name 


N 


M (1O 4 M ) 


Profile 


IMF 


n (pc) 


r v (pc) 


fb 


hcnle5b5 


10 5 


5 


Plummcr 


K91 a [0.1,50] M© 


51 


8.5 


0.05 


hcnle5b!0 


10 5 


5 


Plummcr 


K91 [0.1,50]M Q 


51 


8.5 


0.1 



^Kroupa mass function as described in iKroupa et al.1 (| 199lh . 



-12- 



6 



C\2 

6 



"i 1 1 r 



i 1 1 r 



i 1 1 r 



"i 1 r 




j i i L 



j i i L 



J I I L 



J I L 







5xl0 3 



10 4 
t (Myr) 



1.5xl0 4 



Fig. 3. — Comparison of evolution of the core and the overall binary fraction. Dashed red 
lines show results from lHurley et al.l (120071 . their K100-5 run). Solid black lines show results 
from CMC run using the exact same initial cluster parameters. In both sets the top line show 
the binary fraction within the core, and the bottom line show that for the whole cluster. 
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Fig. 4. — Same as Figure[3j using results from simulation hcnle5bl0. The CMC results are 
compared with the K100-10 simulation of Hurley07. 
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we reproduce the results presented in Hurley07 and verify that the overall /& remains close 
to the primordial value whereas /& )C increases over time. Similar results are found for the 
simulation hcnle5bl0 (Figure HJ). 

We now focus on the evolution of the number of binaries in the core (iV& )C ). The evolution 
of the total number of core binaries is interesting for various reasons. The formation rates of 
interesting stellar objects such as X-ray binaries and blue straggler stars and their properties 
are directly dependent on N hi( . , motivating many detailed studies focusing on the evolution 



ofiV fe.g. lHurlev et al.ll2002l : llvanova et al.ll2005l . 12009 . 12008 : iFregeau et al.ll2009f ). On one 



hand the core binary number (A^ jC ) increases due to mass segregation. On the other hand 
strong interactions involving binary-single (BS) and binary-binary (BB) encounters can lead 
to direct physical collisions or destruction of soft binaries and reduce A^ iC . In addition, 
binary stellar evolution can destroy binaries via evolution-driven mergers and disruptions. 
Since the evolution of Nb tC is dependent on these competing effects it is not simple to predict 
its evolution a-priori. Figure [5] shows the evolution of the number of core binaries (A^ iC ) 
for CMC run hcnle5b5 and direct iV-body run K100-5 from Hurley07 for comparison. The 
evolution of N^ c is reproduced exactly within the numerical fluctuations of the simulations. 
Over time the number of core binaries (iV& >c ) decreases. 

It is also interesting to study the number fraction of binaries and single stars within 
the core compared to the global population. Although A^ >c decreases over time, due to mass 
segregation effects the number of single stars within r c decreases more. Figure M shows the 
evolution of the number fractions of single stars (n S)C ) and binaries (ra?, )C ) within r c for the 
same simulations as above. During the first ~ 10 4 Myr n^^ remains more or less constant 
whereas, n SjC decreases by ~ 0.5 of the initial n SjC due to mass segregation effects. Followed 
by this phase BS/BB interactions as well as stellar evolution destroys core binaries decreasing 
n h c . However, throughout the evolution n h c > n SjC . The combined effects of the above leads 
to the overall increase in ft> tC over time as seen in Figures [3] and HI Note that although 
qualitatively CMC results and the direct iV-body results agree, the agreement is not as 
excellent as the previous comparisons. For example, for the evolution of nt, )C there can be 
upto « 20% difference in the absolute value depending on the age of the simulated cluster. 
The reason behind this larger difference compared to the excellent agreement for the evolution 
of iVb jC (Figure ED originates from the approximations adopted in the tidal treatment in MC 
methods. The criterion based tidal removal of stars adopted in CMC loses stars from the 
tidal boundary at a relatively lower rate (Figure [7]). Hence, at a given time the total number 
of bound single and binary stars in CMC are higher than those in Hurley07 making both rib c 
and n SjC calculated using CMC systematically lower than the same calculated in Hurley07. 

Another interesting result presented in Hurley07 is the evolution of the fraction of bina- 
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Fig. 5. — Evolution of number of binaries within r c (Nj, c ). The solid black line shows the 
results from CMC run hcnle5b5 and the red dashed line shows results from direct iV-body 
run K 100-5 from Hurley07. 
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Fig. 6. — Evolution of the fraction (by numbers) of binary and single stars within the core. 
Thick black lines show the results from CMC. Thin red lines show the same from Hurley07. 
Solid and dashed lines show the number of singles and binaries within r c , respectively. All 
numbers are normalized with the total number of that species at that time in the cluster 
(e.g., Nb^ c /Nb for the binaries). Both results clearly show the effects of mass segregation 
since throughout the evolution a higher fraction of binaries reside in the core. 
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Fig. 7. — Comparison of the evolution of the number of single and binary stars that remain 
bound to the cluster. Thinner lines show results obtained from Hurley07. Thicker lines show 
results from CMC run using the same initial cluster parameters. Solid lines for both cases 
show the number of single stars bound at any given time. Dashed lines show the same for 
binaries. All numbers are normalized to the initial number of the same species (single or 
binary) . 
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Fig. 8. — Evolution of the fraction of double-degenerate binaries within r c with respect to 
all core binaries. Black solid and red dashed lines show results from CMC simulation and 
direct iV-body simulation in Hurley07, respectively. 
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ries in the core where both components are compact objects. We call them double-degenerate 
binaries, following Hurley07. In Figure [8] we show the evolution of the fraction of double- 
degenerate core binaries for CMC run hcnle5b5 and direct iV-body run K100-5. The fraction 
of double-degerate binaries in the core depends on all physical processes in the cluster in a 
complicated way. Two-body relaxation drives mass segregation in the cluster determining 
the densities at different radial regions of the cluster as well as radius dependent velocity 
dispersion. This in turn directly affects the local BS/BB scattering cross-section at a given 
time consequently determining the survivability of a given binary at some radial position in 
the cluster and also the properties of the stellar members in a binary and the binary orbit. 
Changing the binary stellar properties and their orbital properties in turn changes the evolu- 
tionary pathways taken by the binary members and consequently compact object formation. 
Stellar evolution and dynamical effects thus in tandem affect the fraction of double-degerate 
binaries in the core. The excellent agreement between CMC results with Hurley07 convinces 
us that not only the dynamical effects, but also the stellar evolution, and the rate of compact 
object formation, are modeled at least as accurately as in the direct iV-body code NBODY4. 



4. Comparison With Simulations Without Stellar Evolution 

We now examine the effects of stellar evolution on the evolution of the global observable 
properties of a GC, by performing comparisons to simulations without stellar evolution. The 
initial conditions for these simulations are summarized in Table [21 

Our most recent earlier paper, Paper IV, showed results from simulations without stellar 
evolution, but all other physical processes were included. In the absence of an implementation 
of full single and binary stellar evolution Paper IV restricted itself to simulations with a 
narrow range of masses in the IMF. We first compare the results with stellar evolution with 
a small subset of the previous runs from Paper IV without stellar evolution with the narrow 

Table 2. Initial conditions for comparison runs including and leaving out stellar evolution 



Name 


N 


Profile 


IMF 


fb 


kw4b03 


10 5 


King 


KOI [0.1, 1.2] M 


0.03 


kw4bl 


10 5 


King 


K01[0.1,1.2] M 


0.1 


kw4b3 


10 5 


King 


KOI [0.1, 1.2] M 


0.3 


kw7b0 


5 x 10 5 


King 


KOI [0.1, 18.5] M Q 





kw7bl 


5 x 10 5 


King 


KOI [0.1, 18.5] M© 


0.1 
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IMF as an example. Since this is for the purpose of comparison, we use the apocenter 
criterion (§2]) for the tidal treatment to be consistent with Paper IV for these simulations. 

For each of these simulations the initial stellar positions and velocities are chosen from 
a King profile with the concentration parameter Wq = 4. For each sim ulation iV, = 10 5 . The 



stellar IMF is chosen from the stellar MF presented in iKroupal (120011 . henceforth KOI). The 
initial binary fraction is varied between /& = 0.03,0.1, and 0.3. The mass of each binary 
companion is chosen in the range 0.1 — m p from a uniform distribution in mass ratios. The 
binary periods are chosen from a distribution flat in logarithmic a intervals within physical 
limits, where, the hardest binary has a > the sum of the stellar radii of the companions and 
the softest binary is at the local hard-soft boundary. Binary eccentricities are thermal (e.g. 



Heggie &: Hutl 120031 ). For each of these initial conditions one simulation is done including 



stellar evolution and the other leaving it out. 

We find that, even for the simulations with a small range of initial stellar masses, where 
the stellar evolution mass loss is not as severe as in a realistic cluster, for low /& stellar 
evolution can influence the overall cluster evolution to a certain extent. Figure M shows the 
evolution of r c j r^. The results are shown for runs kw4b0 . 03 , 0.1, 0.3 (see Table [2]) . From 
top to bottom the primordial binary fractions /& are 0.03, 0.1, 0.3, respectively. For /& = 0.03 
even with the narrow mass range the two curves start diverging when the most massive stars 
(in this case 1.2 M ) evolve off their MS and lose mass via compact object formation after 
« 3.4 Gyr. 

Binary interactions take place throughout the evolution of the cluster. As the initial /& 
increases, energy available from super-elastic scattering of binaries becomes relatively more 
important compared to the energy produced from stellar evolution mass loss. Thus for this 
narrow range of masses as the binary fraction is increased, the difference between the results 
from simulations including stellar evolution and results without including stellar evolution 
reduces. For example, evolution of the cluster with initial /& = 30% is very similar with and 
without stellar evolution taken into account. The only difference is that at the quasi-steady 
binary-burning phase including stellar evolution makes r c /r^ bigger by about 30%. In each 
of these clusters the central densities are not very high (~ 10 4 M Q /pc 3 ) so direct SS collisions 
are not dominant. When direct SS collisions are more important in a much denser cluster, 



this behavior may change (IChatterjee et al.ll2008l ). 



The difference in the evolution of the global properties depending on whether stellar 
evolution was included or not is, of course, a lot more dramatic when a more realistic IMF 
with a wider mass range is used. Here we use a King profile with central concentration 
parameter Wo = 7. The IMF is according to the K01 stellar MF in the range 0.1 — 18.5 M . 
Two such clusters are simulated one with no primordial binaries and the other with /& = 0.1. 
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4 8 12 

t (Gyr) 

Fig. 9. — Comparison of the evolution of r c /rh with (solid black) and without (dashed red) 
stellar evolution, both using CMC. Results without stellar evolution were already presented 
in Paper IV. Each run starts with 10 5 objects. The velocities and positions of the objects are 
chosen from a King profile with an initial Wo = 4. The masses are chosen from a Salpeter 
MF in the range 0.2 — 1.2 M . From top to bottom the initial binary fractions are 0.03, 0.1, 
and 0.3. 
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Note the dramatic difference in the evolution of the simulated clusters. The evolution 
during the initial ~ 10 2 Myr is dominated by the mass loss via stellar evolution of the high 
mass stars and compact object formation (Figure [TTTl . This phase is clearly distinguished 
by the initial steep expansion of the cluster (Figure [TUT) . This phase is followed by a slow 
contraction phase. In this phase, two-body relaxation drives the evolution. The high mass 
stars have already evolved off the MS and the stars remaining in the cluster are evolving at a 
much slower rate. The transition between the initial stellar evolution driven expansion and 
the slow contraction happens when the energy generation rate from stellar evolution mass 
loss becomes less than the outwards energy diffusion rate from the core due to relaxation. 
The cluster then keeps contracting until the central density increases so much that BS/BB 
interaction rates become high enough and the energy injected by the hard binaries (via 
super-elastic scattering) balances the energy diffusion rate from relaxation. The cluster then 
reaches the binary-burning phase. All these phases are clearly seen in Figure [TD1 bottom 
panel. 

Even with this moderately broad range in mass, the clusters without stellar evolution 
contracts rapidly and are driven towards a quick collapse. If there are primordial binaries, 
the binary-burning phase starts relatively early (~ 1 Gyr, Figure [TUj) . On the other hand, 
when stellar evolution is included, even without any primordial binaries the same cluster 
may still be in the slow contraction phase at Hubble time. For the cluster with primordial 
binaries in this case the binary-burning starts only after 11 Gyr (Figure [TOl) . 



5. Results for Realistic Galactic Globular Clusters 

We have validated CMC by extensive comparisons with direct iV-body results (§3]). 
Moreover, we have shown the importance of including stellar evolution in cluster modeling 
including a realistic stellar IMF (SJ). We now simulate a large grid of clusters with realistic 
initial conditions for 12 Gyr taking all physical processes into account, including primordial 
binaries and single and binary stellar evolution, and the full observed stellar mass range 
spanning three orders of magnitude. Our goal here is to simulate cl usters with realistic 



initia l conditions motivated from observations of young clusters (e.g., IScheepmaker et al. 



20071 ) and find whether at a simulated cluster age of ~ 12 Gyr, a typical age for the GGCs, 
the simulated clusters show similar observable properties (e.g., r c , and r c /r h ) as the observed 
population. 

The proper initial conditions for the GGCs are uncertain, however. Moreover, it is 
hard to infer uniquely the initial conditions from the present day observed cluster prop- 
erties since the observed cluster global properties as well as their galactic orbits can be 
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Fig. 10. — Comparison between two sets of simulations, all using CMC, but in one set stellar 
evolution is included and in the other it is not. Both panels show simulations with an initial 
King profile cluster with Wq = 7 and a Kroupa 2001 IMF in the range 0.1 — 18.5 M Q . The 
top panel shows results with no primordial binaries. The bottom panel shows results with 
0.1 initial /&. For all simulations the initial number of objects is 5 x 10 5 . On both panels 
dashed red lines are for simulations leaving out stellar evolution and solid black lines are for 
simulations including stellar evolution. A dramatic difference is clearly noticeable caused by 
stellar evolution mass loss. 
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Fig. 11. — Rate of stellar evolution mass loss as a function of time. Within 10-100 Myr 
most of the mass is lost due to stellar evolution while high mass stars M* > 5 M evolve off 
their MS. Followed by this initial phase the mass loss rate from winds and compact object 
formation is low. 
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quite uncertain (e.g., iHeggie fc Gierszl 120081 ) . Hence, rather than trying to create a de- 



tailed model for any particular cluster we compare the collective results of all our grid runs 
with the observed GGC properties as a whole . For comp arison the GGC properties are 
extracted from the Harris Catalogue for GGCs (jHarriall996l . and the references therein; also 
see http : / /www.physics.mcmaster.ca/Globular.html[ ). When an observable is not reported 
in the catalogue for a cluster, we exclude that cluster from comparison. In the following 
subsections we explain the initial setup of the grid of simulations and present our results. 



5.1. Initial Conditions 



We simulate clusters with a large grid of initial conditions. All simulated clusters has 
a fixed initial virial radius r v = 4pc (corresponds to an initial pa 3pc). Observations 
indicate that the effective radius of both young and old c l usters are rather insensitive to 



the cluster mass, and metallicity (e.g. lAshman fc Zepfil200ll ; IScheepmaker et al.l 120071 . 120091 ) 
and has a median value of ~ 3pc. In addition, observations of old massive LMC clusters, 
old GCs in NGC 5128, old clusters in M 51, as well as the GGCs indicate that the effective 
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To restrict the huge parameter space to a certain extent we place all our simulated 
clusters in a circular orbit at a Galacto-centric distance of tqc — 8.5 kpc, where the Galactic 
field is not so strong that the tidal stellar loss dominates the cluster's evolution. Choosing 
a circular orbit for the simulated clusters is a simplification, however, the results should 
still be valid for eccentric orbits with some effective pericenter distance of 8.5 kpc (e.g., 
Baumgardt fc Makinoll2003l ). The Galactic tidal field and consequently the initial rt for the 
clusters are calculated using a Galactic rotation speed Vq = 220km/s following the standard 
practice. 

For the set of runs we vary iVj between 4 — 10 x 10 5 , the initial Wo for King models in the 
range 4 — 7.5, and initial /& between — 0.1. For each case we choose the stellar masses of the 
primaries from the MF presented in K01 in the range 0.1 — 100 M©. The masses of each binary 
companion is chosen from a uniform distribution of mass ratios in the range 0.1 — m p M & . 
a for the binaries are chosen from a distribution flat in log within physical limits, namely 
physical contact of the components and the local hard-soft boundary. Although initially each 
binary is hard at its position it may not remain so during the evolution of the cluster. The 
cluster contracts under two-body relaxation and the velocity dispersion increases making 
initially hard binaries soft. Moreover, binaries sink to the core due to mass segregation 
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where the velocity dispersion is higher than the velocity dispersion for the binaries at t — 0. 
We include these soft binaries in our simulations. We let the cluster dynamics disrupt these 
binaries via BS/BB interactions. So at any instant of time soft binaries are allowed in the 
cluster as long as they have not been disrupted naturally via dynamical encounters yet. This 
is closer to reality and this strategy is adopted since soft binarie s can act as an energ y sink 
and can contribute to the overall cluster energetics significantly (IFregeau et al.ll2009l ). 

Each cluster is evolved for 12 Gyr including all physical processes — two body relaxation, 
stellar evolution, strong encounters like BB, BS, and SS collisions. For clusters that reach 
a deep-collapse phase, the CMC time steps become minuscule and the code grinds to a 
halt. We stop our simulations at that point for these clusters. Note that in reality, the 
deep collapse phase is halted via formation of the so called three-body binaries and the 
cluster enters into the gravo-thermal oscillation phase. Since in CMC we do not include the 
possibility of creating new binaries via three body encounters, we do not address this phase 
at this stage. However, this is not a serious limitation for this study since all simulated 
clusters that reach this phase within 12 G yr had a primordial fb = 0, which is not realistic 
(e.g., see most recently iDavis et al.ll2008l ) and simulated for academic interests only. None 
of the simulated clusters enter into the deep-collapse phase before ~ 9 Gyr. The properties 
of all these grid simulations are summarized in Table [31 



5.2. Results 

Here we present some basic observable properties of the simulated clusters and compare 
them with the same properties of the observed GGCs. For each of these comparison plots 
the evolution of a cluster property is shown with the distribution of the same property 
in the GGC population including all GGCs where observation of the concerned property 
exists. Since we restrict the galacto-centric distance of our simulated clusters for this study 
to be 8.5 kpc ( §5.ip we also show the observed distribution of the GGCs with pericenter 
distances from the Galactic center within 7 — 10 kpc to be consistent in the comparisons. 
Note that the purpose for this comparison is simply to ensure that the simulated cluster 
properties agree well with the observed GGC properties. We do not intend to create a 
present day distribution for these properties since for that a probability distribution for the 
initial conditions is required, which is poorly constrained and beyond the scope of this study. 



Note that for the observed GGCs only the sky projection of the characteristic radii such 
as r c and are known. Hence, in order to be consistent in our comparisons with the observed 
population we show the evolution of the 2D projections of r c , and r c /r^ for all simulated 
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Fig. 12. — Evolution of the 2D projection of r c for all simulated clusters. The black, red, 
and blue lines are for clusters with initial /& = 0, 0.05, and 0.1, respectively. A few (6) 
clusters with f b = go into deep-collapse within a Hubble time. We stop integrations for 
those clusters when this phase is reached. The r c values for the observed GGCs are also 
shown in histograms. The solid histogram is for GGCs with pericenter distances between 
7 — lOkpc. When the orbital eccentricities are unknown, a circular orbit is assumed. The 
dashed histogram is for all GGCs where a measurement for r c exist. 
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Fig. 13. — Same as Figure fT2l for the ratio of the 2D projection of the r c to the 2D projection 
of the Th for all simulated clusters. The r c /rh values for the observed GGCs are also shown 
in histograms. The solid and the dashed histograms are from GGC populations selected as 
in Figure [T2"l 
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clusters (e.g., Figures [121 and [13]) . The sky projections for all simulated clusters are done 
assuming spherical symmetry. Both r c as well as the r c /r/j values of the simulated clusters 
agree well with the observed values in the GGC population producing values at 12 Gyr close 
to the peak of the observed distribution. 

We should remind the readers, however, that these r c and values are not exactly the 
quantities observed dire ctly. As mentioned before, r c is the density-weighted core radius 
(jCasertano fc Hutlll985l ). related to a virial radius in the core, co mmonly used in iV-body 
simulations, and can differ from an observed r c by a factor of a few (jHurleyll2007l ). Similarly, 
only the half-light radius is observed which may differ from the half-mass radius of a cluster. 
For example, for a typical simulated cluster elf 3n4 the half light radius including all stars 
is 4.7 pc. If the giant stars are excluded (a common practice for observers) for the same 
cluster the half light radius is 4.1 pc. The theoretically calculated half mass radius for the 
same cluster at the same age is 7pc. 

Nevertheless, one should remember that without including stellar evolution the simu- 
lated r c /r h values including primordial binaries were found t o be about an order of magn itude 



smaller than in the observed population (e.g., Paper IV, IVesperini fc Chernofflll994l ) and 
several studies proposed differen 
large observed rjr h values (e.g 



2008 




Freeeau 


2008 



Mackey et al.l 120081 ) . It is thus quite exciting to find such agreement simply by including 
stellar evolution in the simulations without the need for any fine tuning with the initial 
conditions or exotic scenarios. 

To focus on the distinct evolutionary stages of the clusters we now choose three clusters 
from our large grid of simulations. These clusters are representative of clusters in three 
distinct end stages. Cluster clf3n4 is at the slow contraction phase at the integration 
stopping time and cluster age t c \ = 12 Gyr. Cluster c3f 2nl completes the slow contraction 
phase at t c \ ~ 10 Gyr, reaches the binary-burning quasi-steady phase and remains in the 
binary-burning phase until the integration stopping time t c i = 12 Gyr. Cluster c8flnl 
reaches the deep collapse phase at t cl w 9 Gyr (Figured!]). Integration is stopped after this 
stage is reached. As mentioned earlier, cluster c8f lnl has no primordial binaries and shown 
only as a limiting case for comparison. For each of the simulated clusters the three distinct 
phases of cluster evolution are clearly observed. All simulated clusters first expand due to 
stellar evolution mass loss during the first ~ 1 Gyr. Followed by this initial expansion the 
clusters slowly contracts due to two-body relaxation. This slow contraction phase ends in 
the quasi-steady binary-burning phase for clusters with primordial binaries (Table [3] runs 
except ciflni, where i e [1,4]). Clusters without primordial binaries go into deep-collapse 
directly at the end of slow contraction. 
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Fig. 14. — Evolution of r c /rh (top), r c (middle) and (bottom) for three qualitatively 
different clusters. Results from runs clf3n4, c3f2nl, and c8f lnl are shown in all three 
panels with solid, short-dashed, and long-dashed lines, respectively. Runs clf3n4, c3f 2nl, 
and c8flnl at their final stage of simulation are in the slow contraction, binary-burning, and 
deep-collapse phase, respectively. 
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Fig. 15. — Evolution of the central stellar mass density (p c , top) and the central r.m.s. 
velocity (f c , rms , bottom) for three example runs elf 3n4 (solid), c3f 2nl (short-dashed), and 
c8f lnl (long-dashed). p c and v c ^ Tms decrease sharply during the initial stellar evolution 
driven expansion of the clusters. p c increases slowly during the slow-contraction phase. In 
the binary-burning phase p c attains a quasi-steady value. A sharp increase in p c is observed 
at the deep-collapse phase. v c ^ rms remains more or less steady followed by the initial decrease. 
During deep-collapse v c ^ ms increases sharply. 
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The central density for each simulated cluster first decreases sharply during the initial 
stellar evolution dominated phase due to the early expansion of the core. During the slow 
contraction phase the cluster stellar density increases steadily and reaches a quasi-steady 
value during the binary-burning phase (Figure [TS]). The central velocity dispersion (v c>rms ) 
decreases sharply during the stellar evolution dominated phase. After that v^ rms reaches a 
steady value of ~ 10 km/s. The final value of v^ rrns depends on the evolutionary stage of the 
cluster as well as the total mass in the core. Note that the value of v Ctrms for runs c3f 2nl and 
c8f lnl are similar, since the core masses are comparable, whereas a more massive cluster 
elf 3n4 shows a higher v c ^ rms as expected. v c ^ rms for run c8f lnl starts to diverge from v Cjrms 
for run c3f 2nl only when the former reaches the deep-collapse phase. 

Figure [TBI shows the surface density profiles for the total luminosity and number of stars 
for clusters clf3n4, c3f 2nl, and c8f lnl at the end of simulation. For the first two clusters 
the t c i = 12 Gyr. The third suffers a deep-collapse at ~ 9 Gyr; the profile at the end of 
simulation is shown in that case. We find the best fit single-mass King profile parameters 
minimizing the x 2 statistic from a grid of detailed Kin g models, soly ing the Poisson equation 
where the mass density is calculated self-consistently (]Miocchill2006l . the fitting program was 
kindly provided by Miocchi). Since for old GGCs and similarly for our simulated clusters the 
mass range of the stars at the final stage is narrow, a single-mass King profile is sufficiently 
accurate to predict the cluster parameters such as the King core radius and concentration 
(see Figure [TBI . Furt hermore, we adopt a sing le-mass King fit since observers often follow 



this assumption (e.g.. iDalessandro et al.ll2008l ). The deep-collapsed cluster, c8flnl clearly 



shows a very different projected density profile compared to the other two clusters and cannot 
be represented with a King density profile (Figure [TBI). The collapsed cluster do not have 
a well defined core as seen in the steady increase in the stellar number surface density. For 
cluster elf 3n4, which is in the slow contraction phase, a standard King density profile is an 
excellent representation of the simulated density profile. The density profile in the binary- 
burning cluster c3f 2nl is close to a King profile, however, near the central region there is 
a hint of a power-law density profile expected from observed core-collapsed clusters. In this 
region, a power-law is a better repres entation than a King profile for this cluster ind icating 
a self-similar collapse (Figure [TBI e.g. iHeggie fc Hutll2003l ; iBinney fc Tremaind 120081 ). 



We call the core radius and the concentrations calculated using the best fit King model 
as r c o b s and c Q b s respectively. Table [3] shows a full list of these values for all our simulated 
clusters. However, as shown in Figure [16] these values are not correct for the deep-collapsed 
clusters. Furthermore, for the binary burning clusters a King profile may not be a good fit. 
Nevertheless, most of our simulated clusters are in the slow contraction phase at t d = 12 Gyr, 
where a King density profile is an excellent fit to the data. The luminosity profile is noisy 
due to the presence of a few high luminosity stars (Figure [TBI) . If only stars with a stellar 
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Fig. 16. — Radial profiles for the stellar luminosity density (top) and the stellar number 
density (bottom) at t c \ = 12 Gyr for clusters clf3n4, c3f2nl, and c8flnl, from left to right, 
respectively. The error bars on each panel show the Poisson noise of the data. In the top 
panel the black circles show the luminosity density for all stars taken into account. The red 
squares show the luminosity density taking into account only stars with stellar luminosities 
L* < 20 L . On each panel the dashed blue line shows the best fit King model to the data. 
Bottom left panel is a cluster in slow contraction phase and very well fitted by a King density 
profile. The bottom middle panel is a cluster in the binary burning phase. A King profile 
still works for most parts, however, a hint of a power-law profile is observed (red long dash). 
Bottom right panel is a deep-collapse cluster showing a power-law profile and a King profile 
is a very poor fit. 



-34- 



luminosity < 20 L are taken into account, the profile is less noisy. 

The stellar properties such as the stellar luminosity and radius for each star in the sim- 
ulated cluster are calculated in tandem with the dynamical evolution of the cluster using 
BSE (§2]). From the stellar luminosity and the radius the black-body effective temperature 
can be calculated. Figure [17] shows an example of Hertzprung- Russell diagram (HRD) ob- 
tained from the run clf3n4 at t c \ = 12 Gyr. All binaries for this cluster are assumed to be 
unresolved. The effective temperature of a binary is approximated by a luminosity weighted 
average. Features of a realistic HRD including the MS of the stars, the binary sequence, giant 
branch, and single and binary white dwarf cooling sequences can be clearly seen. Moreover, 
exotic stars such as the BSSs are produced. Here, we define any MSS with a stellar mass 
M* > 1.1 Mto for the cluster as a BSS and find that these stars populate the area of HRD 
expected from observations (Figure [T51) . Here Mto is the MS turn-off mass for the cluster. 
We find that the numbers of BSSs in these clusters depend on the initial conditions as well 
as the evolutionary history. For example, clusters clf3n4, c3f 2nl, and c8f lnl host 52, 16, 
and zero BSSs, respectively, at the time when these snapshots were taken (12 Gyr for the 
first two and integration stopping time ~ 9 Gyr for the third; Figure [TBI) . A more system- 
atic study on the correlations between the total number of BSSs and the cluster observable 
properties is underway. 



6. Summary and Conclusion 

We report the recent update in the development of the Henon based MC code CMC, 
developed at Northweste rn. We have added a fittin g formulae based single and binary stellar 



evolution using BSE by iHurley et al.l (120001 l2002f ) in addition to the already incorporated 



physical processes such as two-body relaxation and strong interaction including BS, BB 
and stellar collisions (Papers I-IV). Thus we are now able to model realistic dense massive 
clusters including all relevant physics with realistic stellar IMFs in our simulations. We 
test the code extensively and compare our results with previously published direct iV-body 
results to validate CMC (DEI). 



In spite of the differences of the basic numerical methods we find that the agreement 
between CMC results and direct iV-body results is excellent (in particular Figures [TJ{8_]) . The 
close reproduction of the evolution of the core fb and the overall /& warrants special mention. 
The evolution of fb is related to all physical processes relevant in the cluster. Two-body 
relaxation drives mass segregation. Binaries being more massive than typical single stars 
mass-segregate towards the center. In the core these binaries interact, and can get destroyed 
via BS/BB interactions. Throughout the evolution the cluster binaries evolve and can merge 




Fig. 17. — An example of a synthetic HRD for run clf3n4 from the set f ixed-7\, at t c i = 
12 Gyr. Each dot is a bound object in the cluster (a single star or a binary). All binaries are 
assumed to be unresolved. The T e ff for a binary is the luminosity weighted temperature. 
The single and binary MSs of the cluster are clearly seen. The giant branch, WD cooling 
sequence, and BSSs are also observed. The stars in between the single MS and the WD 
binary sequence are binaries with MS and WD compact object components. 
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Fig. 18. — Synthetic HRD for the three sample simulated clusters clf3n4, c3f2nl, and 
c8f lnl. For each of the HRDs the region near the tip of the MS is shown. Each dot is 
one bound object (single or binary) in the cluster. The BSSs for each cluster are shown 
as crosses; red and blue crosses denote single and binary BSSs. Here, BSSs are defined as 
the MSSs with mass m* > I.IMto for the cluster at its age. Cluster clf3n4 has a higher 
number of BSSs among the three example clusters (7 singles, 45 binaries) because of this 
cluster's higher iVj and initial /& compared to cluster c3f 2nl (9 singles, 7 binaries). Cluster 
c8f lnl has no BSSs at this age. 
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or disrupt simply through binary evolution. The galactic tidal field tidally strips low mass 
stars from the cluster tidal boundary. Thus obtaining the right evolution of /& indicates all 
these physical processes are implemented correctly. In addition the close agreement between 
the fraction of double-degenerate binaries in the core (Figure [S]) bolsters our belief that not 
only the dynamical processes are modeled correctly, but also stellar evolution and compact 
object formation are at least as accurate as the direct iV-body code NBODY4. 

Although the core properties are obtained accurately using CMC, a larger difference 
is found whenever a quantity involving the total number of bound stars in the cluster is 
compared. For example, the agreement in the evolution of the fractional number of bound 
stars, although is quite good (Figures EJ Ej) given the drastically different methods of simu- 
lations, can differ by at most fa 20%. These differences are dominated by the tidal mass loss 
effects which is hard to model within M C methods and can only be addressed in a criterion 



based way ( §2.3t also see e.g., Paper IV lGiersz fc Spurzeml 120001 ) . A more detailed study in 



characterizing the orbits in a cluster potential and tidal effects is underway, but is beyond 
the scope of this work. 

Our results show that including stellar evolution and a realistic IMF dramatically change 
the evolution of a star cluster and to model a realistic star cluster, inclusion of this 
process is vital. The initial expansion of the cluster driven by stellar evolution mass loss 
significantly prolongs the slow contraction phase. Depending on the initial properties of the 
cluster, even without any primordial binaries the slow contraction phase may last more than 
a Hubble time for clusters typical for the GGCs (e.g., Figures [TU1 IT2"]) . On the contrary, 
exclusion of this important effect leads to a quick contraction of the cluster due to mass 
segregation even if only a moderately broad range of stellar IMF is used (Figure [TO]) . We 
also show that even for simulations with a very narrow stellar mass range, for a relatively low 
fb (Figure ED the inclusion of stellar evolution can give r c /r h values ~ 10% larger compared 
to when it is left out. 

One of the biggest uncertainties of studying the evolution of dense, massive star clusters 
is in determining the initial conditions. The detailed evolution of a cluster depends on 
various initial properties including the initial effective radius, mass, fb, concentration, and 
the galactic tidal field and estimating the initial conditions using the present day observed 
properties of a cluster is hard. Moreover, starting from different initial conditions it is 
possible to achieve very similar present day properties, e.g., r c , and r c /rh (Figures [121 and 
[TBI) . In addition, the observed present day values can also be quite uncertain, especially the 
3D orbit of a cluster in the galactic potential is hard to measure. Furthermore, although 
it may be possible to qualitatively understand individual effects of the various physical 
processes on the observable cluster properties, the collective effect is impossible to judge 
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without actual detailed simulations. Thus to understand a population of dense clusters it is 
required to study a large parameter space and study evolution of these cluster in a realistic 
way including all physics in tandem. 

With the recent improvement to CMC it is now possible to truly scan the full parameter 
space realistically without any loss of generality due to its significantly lower computational 
cost compared to the direct iV-body codes and the accuracy and ability to treat all rele- 
vant physical processes. We have started a detailed study to create a population of realis- 
tic globular clusters, representative of the observed GGCs with a large grid of simulations 
with realistic initial condition s motivated by observations of young massive clusters (e.g. 
Scheepmaker et alJl2007L boogf l. Here we have presented some of these simulations to show 



that rather than creating specific clusters it is now beginning to be possible to create a whole 
population of GGCs using CMC in a star-by-star detail. We show that using observationally 
motivated initial conditions, without any need for fine tuning it is possible to create old 
dense clusters very similar to the observed GGCs (Figures [T2l - [TBI . 



Each star in CMC has realistic stellar properties such as luminosity, radius, and effective 
temperature in addition to the mass and position in the cluster (which is sufficient to follow 
its dynamics). Hence, in addition to the global evolution of the clusters it is possible to 
study individual stellar populations in a cluster. For example, we show synthetic HRDs for 
a few simulated clusters from our grid of simulations. All features, including, e.g., the single 
and binary MS, WD cooling sequence, the giant branch, and BSSs, of a realistic HRD can 
be seen in the synthetic HRD (Figures ITT] [T81 . After this crucial improvement to CMC a 
large array of interesting problems are now accessible. For example, a detailed study of the 
observed GGC BSSs, their properties, and the correlations with various cluster properties is 
underway. 



Table 3. List of simulations: Wq is the central concentration parameter for a King profile f King! [l~966l ) . cluster 



mass 



M is in 10 5 M Q , number of bound cluster objects N is given in 10 5 , central stellar mass density p c is in lO 3 M /pc 3 , r c 
and r/j are in pc, and c is the concentration parameter defined as logio(r c /r t ). All final values are extracted from the 
final snapshot of the simulated clusters. r C)0 6 S and c obs are estimated from a single-mass best fit King model to the 2D 

number density at the final snapshot of the cluster. 
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